Dynamical Linear Response of TDDFT with LDA+ U Functional: 
strongly hybridized Frenkel excitons in NiO 
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Within the framework of time-dependent density- functional theory (TDDFT), we derive the dy- 
namical linear response of LDA+f/ functional and benchmark it on NiO, a prototypical Mott insu- 
lator. Formulated using real-space Wannier functions, our computationally inexpensive framework 
gives detailed insights into the formation of tightly bound Frenkel excitons with reasonable accuracy. 
Specifically, a strong hybridization of multiple excitons is found to significantly modify the exciton 
properties. Furthermore, our study exposes a significant generic limitation of adiabatic approx- 
imation in TDDFT with hybrid functionals and in existing Bethe-Salpeter-equation approaches, 
advocating the necessity of strongly energy-dependent kernels in future development. 

PACS numbers: 71.15.Qe, 78.20.Bh, 71.35.-y, 71.27. +a 
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Computing many-body excitations of weakly to 
strongly interacting materials continues to be a major 
challenge for first-principles studies (l|. Following the 
great achievement of Density Functional Theory (DFT) 
in ground-state calculations using the local density ap- 
proximation (LDA) [2|, ISL time-dependent density func- 
tional theory (TDDFT) 0, [H[ has promised to be an af- 
fordable and accurate theoretical framework to study the 
dynamical linear response of weakly interacting materials 
ranging from finite [(Hsl to extended systems flL lol— -1 1 1 . 
While TDDFT is formally exact in describing the time- 
dependent density, in practice the lack of knowledge of 
the proper functional form of the action functional high- 
lights the immaturity of its full potential in studying dy- 
namics of real materials. For example, there isn't yet an 
affordable first-principles local functional that includes 
proper particle-hole interactions at the two-particle level 
to allow the exciton formation, leavin g th e computation- 
ally expensive perturbation theory [IlJl2| (assuming not 
too strong interactions) or the highly parameterized clus- 
ter model [l3[ the only option. Due to this well-known 
limitation, the applicability of TDDFT (with existing ap- 
proximation) to strongly correlated materials remains ba- 
sically unexplored. 

Recently, a significant improvement in the ground state 
calculations of strongly interacting Mott insulators is 
made via the new LDA+Z7 functional 3, 15], in which 
strong intra-atomic Coulomb interactions are introduced 
at the screened Hartree-Fock level. Similar to the hy- 
brid functionals [16j, the inclusion of non-local exchange 
allows opening of the Mott gap and gives also reason- 
able quasi-particle excitation energies. (Here we do not 
limit TDDFT or DFT to the Kohn-Sham framework, and 
we consider LDA+Z7 and other hybrid functionals as im- 
plicit functionals of density within DFT.) This physically 
motivated approach has led to important understanding 



of orbital and charge ordered systems [17H2Q1- On the 



other hand, to date the study of charge excitation within 
LDA+ U scheme has been very limited and its applica- 
bility to strongly interacting Mott insulators remains un- 
clear. It is thus interesting and timely to explore the dy- 
namical linear response of LDA+ U functional within the 
framework of TDDFT. 

In this Letter, we examine the strength and weakness 
of LDA+ U in describing charge excitations of strongly 
interacting Mott insulators, by developing diagrammati- 
cally the dynamical linear response of LDA+ U functional 
within TDDFT framework (TDLDA+Z7). The resulting 
formula in tackling local excitations connects TDLDA+ U 
(and other hybrid functionals) to a Bethe-Salpeter equa- 
tion (BSE) [2l| with intra-atomic Hartree-Fock kernels. 
The framework is then implemented on the basis of 
symmetry-respecting Wannier functions (WFs) TH, 22], 
which not only dramatically reduce the computation ex- 
pense but also facilitate a comprehensive real-space pic- 
ture of local excitons. The integrated methodologies are 
applied to the study of tightly bound d-d Frenkel exci- 
tons in NiO, a representative Mott insulator. Our di- 
agrammatic approach allows a step-by-step elucidation 
of the effects of different components of the interaction 
kernel. Specifically, multiple tightly bound (by ~6eV) 
excitons form inside the Mott gap and are found to hy- 
bridize strongly to give reasonable excitation energies and 
highly anisotropic q-dependent spectral weights observed 
recently (23[. Furthermore, our results illustrate a serious 
general limitation of adiabatic approximation widely em- 
ployed in most state-of-the-art TD-hybrid methods [lif 
and BSE approaches and advocate the necessity 

of energy dependence in future design of approximations 
within these theoretical frameworks. 

Following the standard procedure [28], the linear re- 
sponse can be derived from the equation of motion 
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FIG. 1: (color online), (a) Dyson equation with LDA+ U 
potential v 3 . (b) Dynamical linear response function x 5 cor- 
relation function L(pifi2]hip2) , and the corresponding BSE 
with LDA+ U kernel. 



(Dyson equation) of Kohn-Sham particle shown diagram- 
matically in Fig. [TJ In addition to LDA potentials, 
site- local (s- local) screened Hartree- Fock interactions and 
the double counting terms among local orbitals are in- 
troduced in LDA+ U. Taking the derivative of LDA+ U 
Green's function with respect to the external potential, 
the response function \-> 



X(xiti;x 2 t 2 ) = 



M;fXL Pl hrMpMt2'Mt2)Mfl h2l 



(1) 

is formulated by creating s- local particle- hole (p-h) pairs 
at position x 2 and time £ 2 , with probability amplitude 
M flh 2 = ^p 2 ( x 2)0/i 2 ( x 2), followed by the propagation 
of the p-h pairs expressed in terms of the p-h correla- 
tion function L, and finally the annihilation of s- local 
p-h pairs at xi and a later time t± with probability am- 
plitude M* x ^ as shown in Fig. QJb). The equation of 
motion of L (BSE) involves irreducible kernel /, includ- 
ing bare Hartree (/h), exchange-correlation / xc , s- local 
screened Hartree-Fock, and double counting terms. The 
s-local Fock interaction directly provides the attraction 
between p-h pairs of orbitals that U is applied to. For 
the study of local d-d excitations, the original /lda and 
the double-counting contributions are meant to counter 
each other approximately and should be neglected, leav- 
ing only the s-local screened Hartree-Fock contribution. 
(Besides, the long-range screening of /h is inefficient at 
short distance.) This diagrammatic representation of the 
equation of motion makes it apparent that one is allowed 
to visualize the physical effects of interaction kernels by 
turning them on one after another, as presented below. 

For the study of local excitations, it is most convenient 
and efficient to employ a real-space s-local basis [29] . Fur- 
thermore, it is advantageous to have orbitals that diag- 
onalize the one-particle density matrix such that they 



enter only as either pure particle or hole orbitals. To this 
end, we implemented the above TDLDA+ U framework 
using the energy-resolved symmetry-respecting Wannier 
functions [ljj [22[ that are constructed without mixing 
the occupied and unoccupied bands [30]. In this ba- 
sis, the complicated six-dimensional exciton wave func- 
tion, ^(xx ; ), can be decomposed into small number 
of "bare" exciton wave functions as direct products of 
one particle (<p p ) and one hole (<frh) orbitals, ^(xx') = 
^2ph c ph^p(' x ) ( t ) h( yi ')' This gives an easy visualization of 
exciton wave functions and a direct computation of var- 
ious experimental "form factor". Moreover, the strong 
binding suppresses significantly the kinetics of the exci- 
ton such that only very few short-range neighbors are 
necessary in solving the BSE, significantly reducing the 
computational expense. In fact, even a single-site calcu- 
lation is found quite accurate for NiO [30| . 

The LDA+ U calculation for type-II antiferromagnetic 
(AFM) NiO ({7=8eV, J=0.95eV) is performed by an 
all-electron full-potential method using linearized aug- 
mented plane wave basis [3l|. It gives the correct 
ground state of NiO with high-spin configuration, leav- 
ing only the spin-minority channel active for s-local d-d 
charge excitations. The rhombohedral symmetry, dic- 
tated by the AFM order along [111] direction, splits the 
d- Wannier orbitals into two unoccupied e g , two occupied 
e' g , and one occupied a g orbitals. In the rest of the Let- 
ter, we denote l->(e g i,a g ), V ^(e g2 ,a g ), 2->(e g i,e' gl ), 
2' -+(e 9 2,e g2 ), 3^(e g2 ,e gl ), and 3' -+(e gl ,e' g2 ) for dif- 
ferent p-h pairs (p,h) forming bare excitons. Within this 
notation, xCq?^) 1S thus . M*^L^-(q, uS)M^. Here 
= J e*^' x M x <ix gives the probability amplitudes in 
momentum space corresponding to the anisotro pic form 
factor of inelastic X-ray and electron scattering |23l l32| . 

We now demonstrate in Fig. [2] the physical effects of 
the above interaction kernels on forming the low-energy 
Frenkel excitons in NiO, by switching them on one af- 
ter another in the corresponding BSE of TDLDA+ U. (A 
broadening of r]=0. 2eV is introduced in calculating L in 
accordance to the experimental resolution 0.) As a ref- 
erence, the unbound p-h excitations, Lq, containing no 
interaction kernel, is shown in Fig. [SJa). Lq gives ba- 
sically the creations of p-h pairs across sites (inter-site 
d-d excitations) that contribute at small momentum but 
are eventually overwhelmed by the on-site excitations at 
larger momentum. Not surprisingly, the excitation en- 
ergies correspond to the energy difference of LDA+/7 
eigenvalues encapsulating a large Coulomb interaction 
(of order U) due to the addition of a particle and a hole 
at different atoms. In addition, a broad spectrum (5-10 
eV) is obtained, reflecting the high mobility of such un- 
bound p-h pairs and their strong coupling to the oxygen 
p-orbitals. 

For intra-site excitations, however, the interaction ker- 
nel can dramatically modify the excitations. As shown 
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in Fig. [2(a), upon switching on the Fock kernel within 
the same pairs, which encapsulates the strong binding of 
p-h pairs, the excitations collapse into three well-defined 
doubly-degenerate bare excitons at l-2eV via binding en- 
ergies ~6eV! This seemingly large binding energy is a 
necessary consequence of the large s-local Coulomb re- 
pulsion present in NiO and other strongly correlated ma- 
terials. Preserving the particle number at each site, the 
s-local charge excitations should not be subject to the 
Hubbard U energy scale of adding or removing one par- 
ticle. On the other hand, since one-particle Green's func- 
tions used to build Lo includes the energy scale of [/, such 
effects must be countered by a strong attraction of the 
same scale in any proper theory of p-h pairs, as demon- 
strated here. Benefiting from the lack of decay process 
deep inside the Mott gap, these bare excitons have very 
long lifetime (negligible line width beyond the experi- 
mental resolution). Due to the well-defined point-group 
symmetry of these Wannier orbitals (and thus M^), the 
resulting bare excitons possess well-defined angular struc- 
ture in their spectral weight. As shown in Fig. [2(d), a 
prominent anisotropy of x 1S found with the strongest 
weight along [111] directions and vanishing weight along 
[001] directions, in agreement with the recent experi- 
ment [23j| . In addition, the dipole- forbidden nature is 
revealed clearly by a hollow center at q — >• 0. 

Next, we activate the scattering process among differ- 
ent p-h pairs in the remaining part of the Fock kernel, 
I F . These additional couplings between bare excitons in- 
troduce an interesting effect of strong hybridization and 
split the exciton energies into two sets. As illustrated in 
Fig. [2(b) , this strong hybridization results in new large 
off-diagonal elements of L in the bare exciton basis, and 
the double-peak structure in the diagonal elements. From 
the conventional eigenvector point of view, this is equiva- 
lent to having exciton eigen-functions of each energy con- 
taining superposition of the above bare excitons. Indeed, 
the q-dependence of excitons for the low-energy and the 
high-energy peaks now shows significant differences, as 
shown in Fig. [2(e) and (f). Along [111] directions, the 
former has a vanishing weight (new nodal directions), 
while the latter shows strong enhancement, reflecting the 
anti-bonding and bonding nature of the superposition, 
respectively. Such strong hybridization of multi-excitons 
is absent in weakly correlated systems, but should be ex- 
pected from most strongly correlated systems with open 
shells. 

Finally, upon addition of the s-local Hartree kernel, 
I H , the overall exciton spectrum is further modified. As 
shown in Fig. [2(c), an overall blue shift of the spectrum 
is observed, originated from the screening (weakening) 
of the exciton binding energy via the diagonal elements 
of I H . In addition, without altering the q-dependent 
spectral weight significantly, the off-diagonal elements of 
I H further split the exciton energies. Both features are 
expected from the Kramers-Kronig relation in the typical 




4 

3 

> 2 
5, 1 


-1 



. 1 1 1 1 II 1 1 

-(b) A A 


1 1 1 1 1 1 1 1 
i f - 




\ + h i Pi - 




\ h 2 p 2 - 


~. 1 1 1 1 1 1 1 1 


, , , , 1 , , , , 1 , . 7 



1 1 1 1 1 1 1 1 1 1 

-(c) A 


1 1 1 1 1 1 1 1 1 1 1 1 
A IH - 






1 \\ + ': 




h 2 ° p 2 _ 


~ . . . 1 , , y , 1 


M 1 1 1 1 1 1 1 1 M 



2 

co(eV) 




FIG. 2: (color online). Imaginary parts of p-h correlation 
function Lij denoted as (i;j) by solving BSE for (a) Lo: un- 
bound p-h pair, and L: bare exciton, (b) L with s-local Fock, 
and (c) L with s-local Hartree-Fock. Experimental excitation 
energies are indicated by the arrows, (d), (e), and (f) illus- 
trate the spectral weight of Lo, low-energy exciton in L, and 
high-energy exciton in L, respectively, via 3D isovalue con- 
tours in q-space. The black solid circles indicate q=7A~ 1 . 
Spectral weight in (e)/(f) resembles very well the experimen- 
tal ones corresponding to the blue/red arrows. 



dynamical screening process. 

So far, the above analysis has been performed with 
Wannier functions associated with only one Ni site, which 
makes the proposed framework extremely inexpensive 
and allows a very clean real-space physical picture. More 
extensive calculations [30] including non-s-local (beyond 
one Ni site) effects resulted in a negligible energy reduc- 
tion (<0.1eV) of the overall spectrum, and, in perfect 
agreement with the experiment [23|, a negligible disper- 
sion. Indeed, the Frenkel excitons in NiO are so tightly 
bound and so heavy to propagate that the physics can 
be captured quite completely even with only a single site. 
This further advocates strongly our choice of employing 
Wannier basis for strongly correlated materials in gen- 

3 EES, 0. 



eral 

Overall, compared with the recently measured local 
d-d excitations [23|, [32|, the resulting theoretical spec- 
trum of above TDLDA+ U performs unexpectedly well 
for such a simple approximation. In great contrast to 
the qualitative failure of TDLDA ^ or the RPA re- 
sponse of LDA+£7 that give not even a hint of ex- 
citons, TDLDA+£7 generates successfully the tightly 
bound Frenkel excitons deep inside the Mott gap with 
good energy scales at 1-3 eV. Considering the enormous 
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binding energies of ~6eV and the complexity of the multi- 
exciton coupling, this degree of agreement is quite im- 
pressive, not to mention the good q-dependent spectral 
weight that reflects the high quality of the exciton wave 
functions. Since the analytical structure of LDA+ U is 
the same as the hybrid functionals and its equation of 
motion contains similar kernel to that of existing GW- 



BSE [25|, [26| approximation, the same degree of success 
should be expected from these approaches as well, so 
should above general physical picture of the formation 
of Frenkel excitons. 

On the other hand, our study also reveals a generic 
limitation of TDLDA+ U and all these state-of-the-art 
approximations. While TDLDA+£7 puts the excitons in 
the right energy range, it is unable to reproduce the fine 
structures of the experimental excitons satisfactorily. In- 
deed, while three charge excitation energies in the Mott 
gap were found experimentally (c.f.: Fig. 2), TDLDA+Z7 
only manages to produce two distinct energies. It turns 
out that the inability to split resulting peaks into the fine 
structures is dictated by the analytical structure of the 
approximations. This can be elucidated by reformulating 
rigorously the BSE into an effective Hamiltonian via the 
creation operator of bare excitons b\\ 



i i^j 



where q denotes the energy of the bare excitons in 
Fig.[2fa). In this representation, all the above physical ef- 
fects are made transparent: the hybridization due to the 
off-diagonal elements of I F and I H , and the blue shift 
of the excitonic energy via the diagonal elements of I H . 
This representation also proves that without terms that 
encapsulate inter-boson interactions, e.g. if^, j,b\bjbj'bi> , 
this class of approximations is only capable of produc- 
ing excitons in number equal to that of the bare excitons 
(3-2 = 6 for NiO), subject to degeneracy under the point- 
group symmetry. Thus, this generic limitation can only 
be removed via higher order kernel in BSE or new gener- 
ations of hybrid-functional kernels that explicitly include 
time-dependence. 

This general conclusion is actually more stringent than 
those in the present TDDFT literature. Without uj- 
dependent J, our results demonstrate exciton formation 
with correct scale of large binding energy. Clearly, in this 
crudest level of approximation, the known necessity of 
non-adiabatic (memory-dependent) kernel / xc of TDDFT 



in the Kohn-Sham framework 31 36] would originate 
solely from the spatial reduction of non-local J xc into 
/xc [121 • On the other hand, to allow any fine ("multi- 
plets") structure in strongly interacting systems, a non- 
adiabatic kernel is absolutely necessary even without the 
spatial reduction. Obviously, this is one key aspect that 
almost all the existing approximate functionals lack and 
presents an essential and necessary step toward a proper 



description of local excitations in strongly interacting sys- 
tems, within all the existing theoretical frameworks. 

In summary, a TDLDA+ U method is derived within 
TDDFT framework and implemented on the basis of 
Wannier functions. As a benchmark and an illustration, 
the formation of tightly bound excitons in NiO are ana- 
lyzed step-by-step. A strong hybridization of multiple ex- 
citons is found essential to give reasonable excitation en- 
ergies with a large ~6eV binding energy and the observed 
highly anisotropic spectral weight. Our computationally 
inexpensive approach not only provides detailed insights 
into the formation of Frenkel excitons, but also demon- 
strate an intrinsic non-adiabaticity of the interaction ker- 
nel in allowing fine-structures in the excitation spectra of 
strongly interacting systems. The lack of such essential 
feature in nowadays hybrid- functional and BSE methods 
advocates the inclusion of explicit time-dependence in fu- 
ture design of approximate functionals. 
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ROC for support, and also NCHC of ROC for CPU time. 



[1] 
[2] 

[3] 
[4] 

[5] 
[6] 

m 

[8] 
[9 

[io; 
[ii 

[12 

[13 

[14 
[15 
[16 
[17 
[18 
[19 
[20' 
[21 

[22 
[23 
[24 
[25 
[26 



G. Onida et al, Rev. Mod. Phys. 74, 601 (2002). 

P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 
(1964). 

W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965). 

E. Runge and E. K. U. Gross, Phys. Rev. Lett. 52, 997 
(1984). 

R. van Leeuwen, Phys. Rev. Lett. 82, 3863 (1999). 
M. Petersilka et al., Phys. Rev. Lett. 76, 1212 (1996). 
I. Vasiliev et al., Phys. Rev. Lett. 82, 1919 (1999). 
K. Yabana and G. F. Bertsch, Int. J. Quantum Chem. 
75, 55 (1999). 

N. E. Maddocks et al., Europhys. Lett. 27, 681 (1994). 
W. Ku and A. G. Eguiluz, Phys. Rev. Lett. 82, 2350 
(1999). 

W. Ku et al., Phys. Rev. Lett. 88, 057001 (2002). 

I. V. Tokatly and O. Pankratov, Phys. Rev. Lett. 86, 

2078 (2001). 

M. W. Haverkort et al., Phys. Rev. Lett. 99, 257401 
(2007). 

V. I. Anisimov et al., Phys. Rev. B 44, 943 (1991). 

A. B. Shick et al., Phys. Rev. B 60, 10763 (1999). 

F. Tran et al., Phys. Rev. B 74, 155108 (2006). 

I. Leonov et al., Phys. Rev. Lett. 93, 146404 (2004). 

H. -T. Jeng et al., Phys. Rev. Lett. 93, 156403 (2004). 
W.-G. Yin et al., Phys. Rev. Lett. 96, 116405 (2006). 

D. Volja et al., Europhys. Lett. 89, 27008 (2010). 

E. E. Salpeter and H. A. Bethe, Phys. Rev. 84, 1232 
(1951). 

W. Ku et al., Phys. Rev. Lett. 89, 167204 (2002). 

B. C. Larson et al., Phys. Rev. Lett. 99, 026401 (2007). 
J. Paier et al., Phys. Rev. B 78, 121201(R) (2008). 

G. Onida et al., Phys. Rev. Lett. 75, 818 (1995). 

M. Rohlfing and S. G. Louie, Phys. Rev. Lett. 81, 2312 



5 



(1998). 

[27] L. X. Benedict et al., Phys. Rev. Lett. 80, 4514 (1998). 
[28] G. Baym and L. P. Kadanoff, Phys. Rev. 124, 287 (1961). 
[29] J. A. Soininen et al., Phys. Rev. B 72, 045136 (2005). 
[30] See EPAPS Document No. for supplementary material. 
[31] K. Schwarz et al., Comput. Phys. Commun. 147, 71 
(2002). 



[32] F. Miiller and S. Hiifner, Phys. Rev. B 78, 085438 (2008). 
[33] V. Turkowski et al., Phys. Rev. B 79, 233201 (2009). 
[34] Y. Kim and A. Gorling, Phys. Rev. B 66, 35114 (2002). 
[35] M. Hellgren and U. von Barth, Phys. Rev. B 78, 115107 
(2008). 

[36] N. T. Maitra et al., Phys. Rev. Lett. 89, 023002 (2002). 



